/**
 * @file   Mesh.c
 * @author xlhu <xlhu@Nano>
 * @date   Tue Jan 15 23:21:11 2013
 * 
 * @brief  Various Mesh IO functions
 * 
 * 
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "Mesh.h"
#include "SparseMatrix.h"
#include "constants.h"

/* // the template geometrical variabled defined in constants.c */
/* extern int Edge_end[6][2];     // local edge-vertex index on the reference tet */
/* extern int Face_end[4][3];     // local face-vertex index on the reference tet */
/* extern int Face_edge[4][3];    // local face-edge   index on the reference tet */

/* extern int BM_OUTER; */
/* extern int BM_INTERFACE; */

/* extern int Edge_end[6][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}}; */

/* extern int Face_end[4][3] = {{0,1,2},{0,1,3},{0,2,3},{1,2,3}}; */

/* extern int Face_edge[4][3] = {{0,1,3},{0,2,4},{1,2,5},{3,4,5}}; */

/* extern int Face_edge_end[3][2] = {{0,1}, {0,2}, {1,2}}; */

/* extern int BM_OUTER = 100; */
/* extern int BM_MPI   = 200; */
/* extern int BM_GHOST = 300; */
/* extern int BM_INTERFACE   = 10; */
/* extern int BM_INTERFACE_A = 20; */
/* extern int BM_INTERFACE_B = 30; */



/////////////////////////////////
///  read the mesh head in Original Format
/////////////////////////////////
void readMeshHeadA(char *meshPath,
		   int *N_node, int *N_edge, int *N_element,
		   int *N_point, int *N_line, int *N_triangle)
{
  char fileName[60];
  FILE *fout1;

  sprintf(fileName, "%s/N_node.txt", meshPath);
  fout1=fopen(fileName,"r");
  fscanf(fout1,"%d",N_node);
  fclose(fout1);
  
  sprintf(fileName,""); // empty the string
  sprintf(fileName, "%s/N_edge.txt", meshPath);
  fout1=fopen(fileName,"r");
  fscanf(fout1,"%d",N_edge);
  fclose(fout1);
  
  sprintf(fileName,""); // empty the string
  sprintf(fileName, "%s/N_element.txt", meshPath);
  fout1=fopen(fileName,"r");
  fscanf(fout1,"%d",N_element);
  fclose(fout1);

  sprintf(fileName,""); // empty the string
  sprintf(fileName, "%s/N_point.txt", meshPath);
  fout1=fopen(fileName,"r");
  fscanf(fout1,"%d",N_point);
  fclose(fout1);
  
  sprintf(fileName,""); // empty the string
  sprintf(fileName, "%s/N_line.txt", meshPath);
  fout1=fopen(fileName, "r");
  fscanf(fout1,"%d",N_line);
  fclose(fout1); 
  
  sprintf(fileName,""); // empty the string
  sprintf(fileName, "%s/N_triangle.txt", meshPath);
  fout1=fopen(fileName, "r");
  fscanf(fout1,"%d",N_triangle);
  fclose(fout1);

}

/**
 * Read the original mesh data
 * 
 */

void readMeshDataA(char  *meshPath,
		  double *Node,     int N_node,     //
		  int    *Edge,     int N_edge,     //
		  int    *Element,  int N_element,  //
		  int    *Point,    int N_point,    //
		  int    *Line,     int N_line,
		  int    *Triangle, int N_triangle,
		  int    *Triangle_Element,   // for find neighbour
		  int    *Element_Edge,        // edge-element 
		  int    *Triangle_Edge,       // boundary triangle
		  int    *Line_Edge,
		  int     verbose)           // boundary dof
{
  char fileName[60];
  FILE *fout1;
  int i,idx,dummy;
  double x,y,z;
  int j,k,l,m,n;
  
  if (verbose) printf("Loading Meshes ... \n");   
  // 3. Read the real mesh data
  sprintf(fileName,""); // empty the string
  sprintf(fileName, "%s/Node.txt", meshPath);
  fout1=fopen(fileName,"r");
  for(i = 0, idx = 0; i < N_node; i++, idx += 3){
    fscanf(fout1, "%lf %lf %lf", &x, &y, &z);
    Node[idx+0] = x;
    Node[idx+1] = y;
    Node[idx+2] = z;
  } 
  fclose(fout1); 

  sprintf(fileName,""); // empty the string
  sprintf(fileName, "%s/Edge.txt", meshPath);
  fout1=fopen(fileName,"r");
  for(i = 0, idx = 0; i < N_edge; i++, idx += 2){
    fscanf(fout1,"%d %d", &m, &n);
    Edge[idx+0] = m - 1;   ///  in memory, the id
    Edge[idx+1] = n - 1;   ///  begin from 0!!
  }
  fclose(fout1); 


  sprintf(fileName,""); // empty the string
  sprintf(fileName, "%s/Element.txt", meshPath);
  fout1=fopen(fileName,"r");
  for(i = 0, idx = 0; i < N_element; i++, idx += 4){
    fscanf(fout1,"%d %d %d %d", &j, &k, &m, &n);
    Element[idx+0] = j - 1;
    Element[idx+1] = k - 1;
    Element[idx+2] = m - 1;
    Element[idx+3] = n - 1;
  } 
  fclose(fout1); 


  sprintf(fileName,""); // empty the string
  sprintf(fileName, "%s/Point.txt", meshPath);
  fout1=fopen(fileName,"r");
  for(i = 0; i < N_point; i++){ 
     fscanf(fout1,"%d",&l);
     Point[i] = l - 1;
  } 
  fclose(fout1); 

  sprintf(fileName,""); // empty the string
  sprintf(fileName, "%s/Line.txt", meshPath);
  fout1=fopen(fileName,"r");
  for(i = 0, idx = 0; i < N_line; i++, idx += 2){
    fscanf(fout1,"%d %d", &j, &k);
    Line[idx+0] = j - 1;
    Line[idx+1] = k - 1;
  } 
  fclose(fout1); 

  sprintf(fileName,""); // empty the string
  sprintf(fileName, "%s/Triangle.txt", meshPath);
  fout1=fopen(fileName,"r");
  for(i = 0, idx = 0; i < N_triangle; i++, idx += 3){
    fscanf(fout1,"%d %d %d", &j, &k, &l);
    Triangle[idx+0] = j - 1;
    Triangle[idx+1] = k - 1;
    Triangle[idx+2] = l - 1;
  } 
  fclose(fout1); 

  sprintf(fileName,""); // empty the string
  sprintf(fileName, "%s/Triangle_Element.txt", meshPath);
  fout1=fopen(fileName,"r");
  for(i = 0, idx = 0; i < N_triangle; i++, idx += 2){
    fscanf(fout1,"%d %d", &j, &k);
    Triangle_Element[idx+0] = j - 1;
    Triangle_Element[idx+1] = k - 1;
  } 
  fclose(fout1);


  sprintf(fileName,""); // empty the string
  sprintf(fileName, "%s/Line_Edge.txt", meshPath);
  fout1=fopen(fileName,"r");
  for(i = 0; i < N_line; i++){ 
     fscanf(fout1,"%d",&l);
     Line_Edge[i] = l - 1;
  }
  fclose(fout1);

  sprintf(fileName,""); // empty the string
  sprintf(fileName, "%s/Triangle_Edge.txt", meshPath);
  fout1=fopen(fileName,"r");
  for(i = 0, idx = 0; i < N_triangle; i++, idx += 3){
    fscanf(fout1,"%d %d %d", &j, &k, &l);
    Triangle_Edge[idx+0] = j - 1;
    Triangle_Edge[idx+1] = k - 1;
    Triangle_Edge[idx+2] = l - 1;
  }  
  fclose(fout1); 

  sprintf(fileName,""); // empty the string
  sprintf(fileName, "%s/Element_Edge.txt", meshPath);
  fout1=fopen(fileName,"r");
  for(i = 0, idx = 0; i < N_element; i++, idx += 6){
    fscanf(fout1,"%d %d %d %d %d %d", &j, &k, &l, &m, &n, &dummy);
    Element_Edge[idx+0] = j - 1;
    Element_Edge[idx+1] = k - 1;
    Element_Edge[idx+2] = l - 1;
    Element_Edge[idx+3] = m - 1;
    Element_Edge[idx+4] = n - 1;
    Element_Edge[idx+5] = dummy - 1; // ugly, you can modify it!
  } 
  fclose(fout1);  

  if(verbose) 
    printf("Read %d Nodes, %d Edges and %d Elements in directory %s in the internal format!\n", N_node, N_edge, N_element, meshPath);

}


/**
 * Read Internal mesh format head
 * 
 */
void readMeshHead(char *meshFileName,
		  int  *N_node,
		  int  *N_edge,
		  int  *N_triangle,
		  int  *N_element,
		  int   verbose)
{
  char meshFile[60], dummy[60];
  FILE *f_in;
  sprintf(meshFile,"%s.mesh",meshFileName);
  f_in = fopen(meshFile, "r");
  //skip the first line
  fscanf(f_in, "%s", dummy);
  fscanf(f_in, "%d %d %d %d", N_node, N_edge, N_triangle, N_element); 
  fclose(f_in);

  if(verbose)
    printf("There are %d nods, %d edges, %d boundary tiangles and %d elements.\n", 
	   *N_node, *N_edge, *N_triangle, *N_element);

}

/**
 * Read the true mesh data
 * 
 */

void readMeshData(char *meshFileName,
		  double *Node,    int N_node, 
		  int    *Edge,    int N_edge,
		  int    *Element,        // node dof for the element
		  int    *ElementEdge,    // edge dof for the element
		  int    *ElementMark,    // ElementMark for the element
		  int     N_element,  
		  Face   *TriangleList, int N_triangle,   
		  int verbos)           // boundary dof
{
  char meshFile[60], dummy[60];
  FILE *f_in;
  int i, j,k,l, m,n, p,q,r, s,t, idx, nn;

  sprintf(meshFile,"%s.mesh",meshFileName);
  if(verbos)
    printf("Readding the mesh data in internal format from %s ...\n", meshFile);

  f_in = fopen(meshFile, "r");

  // 0. write the mesh file head, N_triangle is the boundary triangles only

  fscanf(f_in, "%s", dummy);
  fscanf(f_in, "%d %d %d %d", &m, &n, &p, &q); 
  fscanf(f_in, "%s", dummy);
  fscanf(f_in, "%s", dummy);

  // 1. read the nodes' coordinates
  fscanf(f_in, "%d", &nn);
  if(N_node != nn){ 
    printf("N_node is not equation to the file.\n"); 
    return;
  }
  for(i = 0; i < N_node; i++){
    fscanf(f_in, "%lf %lf %lf", 
	   &Node[3*i+0], &Node[3*i+1], &Node[3*i+2]);
  }

  // 2. read 1d element: the edges
  fscanf(f_in, "%d", &nn);
  if(N_edge != nn){ 
    printf("N_edge is not equation to the file.\n"); 
    return;
  }
  for(i = 0; i < N_edge; i++){
    fscanf(f_in, "%d %d", &m,&n);
    Edge[2*i+0] = m - 1;
    Edge[2*i+1] = n - 1;
  }

  // 3. read 2d element: the triangle
  fscanf(f_in, "%d", &nn);
  if(N_triangle != nn){ 
    printf("N_triangle is not equation to the file.\n"); 
    return;
  }
  Face *the_face;
  for(i = 0; i < N_triangle; i++){
    fscanf(f_in, "%d  %d %d %d  %d %d %d  %d %d  %d",
	   &nn,  &j, &k, &l,  &p, &q, &r,  &m, &n, &idx); 

    the_face = (Face *) malloc(sizeof(Face));
    the_face->_id    = nn - 1;
    the_face->_vtx[0] = j - 1;
    the_face->_vtx[1] = k - 1;
    the_face->_vtx[2] = l - 1;

    the_face->_bnd[0] = p - 1;
    the_face->_bnd[1] = q - 1;
    the_face->_bnd[2] = r - 1;

    the_face->_nei[0] = m - 1;
    the_face->_nei[1] = n - 1;

    the_face->_boundary_mark = idx;

    the_face->_next = TriangleList->_next;
    TriangleList->_next = the_face;
  }

  // 4. write 3d element: the tetrahedron
  fscanf(f_in, "%d", &nn);
  if(N_element != nn){ 
    printf("N_element is not equation to the file.\n"); 
    return;
  }
  for(nn = 0; nn < N_element; nn++){
    fscanf(f_in, "%d  %d %d %d %d  %d %d %d %d %d %d  %d", 
	   &idx, &i, &j, &k, &l, &m, &n, &p, &q, &s, &t, &r);
 
    Element[4*nn+0] = i - 1;
    Element[4*nn+1] = j - 1;     
    Element[4*nn+2] = k - 1;
    Element[4*nn+3] = l - 1;
    
    ElementEdge[6*nn+0] = m - 1;
    ElementEdge[6*nn+1] = n - 1;
    ElementEdge[6*nn+2] = p - 1;
    ElementEdge[6*nn+3] = q - 1;
    ElementEdge[6*nn+4] = s - 1;
    ElementEdge[6*nn+5] = t - 1;

    ElementMark[nn] = r;
  }

  if(verbos)
    printf("Read %d Nodes, %d Edges, %d Tets and %d boundary Tris from %s.\n", 
	   N_node, N_edge, N_element, N_triangle, meshFile);

  // 3.end of the output
  fclose(f_in);
  
}

/**
 * Write down all the necessary mesh data information
 * int the Internal Format
 */
void writeMesh(char   *meshFileName,
	       double *Node,    int N_node, 
	       int    *Edge,    int N_edge,
	       int    *Element,        // node dof for the element
	       int    *ElementEdge,    // edge dof for the element
	       int    *ElementMark,    // ElementMark for the element
	       int     N_element,  
	       Face   TriangleList, int N_triangle,   
	       int verbos)           // boundary dof
{
  char meshFile[60];
  FILE *fout;
  int i,j, idx;
  double x,y,z;

  sprintf(meshFile,"%s.mesh",meshFileName);
  if(verbos)
    printf("Recording the mesh data in internal format into %s...\n", meshFile);

  fout = fopen(meshFile, "w");

  // 0. write the mesh file head, N_triangle is the boundary triangles only

  fprintf(fout, "<Head>\n");
  fprintf(fout, "%d %d %d %d\n", N_node, N_edge, N_triangle, N_element); 
  fprintf(fout, "</Head>\n");

  fprintf(fout, "\n");  // separate line
  fprintf(fout, "<Body>\n");

  if(Node == NULL || Element == NULL){
    printf("The Mesh is empty!\n");
    fprintf(fout, "</Body>\n");
    return;
  }

  // 1. write the nodes' coordinates
  fprintf(fout, "%d\n", N_node);
  for(i = 0; i < N_node; i++){
    idx = 3*i;
    fprintf(fout, "%lf %lf %lf \n", 
	    Node[idx], Node[idx+1], Node[idx+2]);
  }
  fprintf(fout, "\n");   // empty line


  // 2. write 1d element: the edges
  fprintf(fout, "%d\n", N_edge);
  for(i = 0; i < N_edge; i++){
    fprintf(fout, "%d %d \n", 
	    Edge[2*i+0] + 1, Edge[2*i+1] + 1);
  }
  fprintf(fout, "\n");   // empty line

  // 3. write 2d element: the triangle
  fprintf(fout, "%d\n", N_triangle);
  int count = 0;
  Face *the_face = TriangleList._next;   // the head node is empty!
  while(the_face != NULL){
    fprintf(fout, "%d \t%d %d %d \t%d %d %d \t %d %d \t%d\n",
	    the_face->_id+1, 
	    the_face->_vtx[0]+1, the_face->_vtx[1]+1, the_face->_vtx[2]+1,
	    the_face->_bnd[0]+1, the_face->_bnd[1]+1, the_face->_bnd[2]+1, 
	    the_face->_nei[0]+1, the_face->_nei[1]+1,  // be 0 if boundary triangle
	    the_face->_boundary_mark);
    the_face = the_face->_next;  // go to the next triangle
    count ++;
  }
  fprintf(fout, "\n");   // empty line

  // check the validaty
  if(count != N_triangle) 
    printf("Error in writing the triangle data into internal format!\n");

  // 4. write 3d element: the tetrahedron
  fprintf(fout, "%d\n", N_element);
  for(i = 0; i < N_element; i++){
    fprintf(fout, "%d \t%d %d %d %d \t%d %d %d %d %d %d \t%d\n",
	    i+1,
	    Element[4*i+0] + 1,     Element[4*i+1] + 1,     
	    Element[4*i+2] + 1,     Element[4*i+3] + 1,
	    ElementEdge[6*i+0] + 1, ElementEdge[6*i+1] + 1, 
	    ElementEdge[6*i+2] + 1, ElementEdge[6*i+3] + 1,
	    ElementEdge[6*i+4] + 1, ElementEdge[6*i+5] + 1,
	    ElementMark[i]);
  }
  fprintf(fout, "\n");   // empty line

  fprintf(fout, "</Body>\n");

  if(verbos)
    printf("%d Nodes, %d Edges, %d Tets and %d boundary Tris written into %s!\n", 
	   N_node, N_edge, N_element, N_triangle, meshFileName);

  // 3.end of the output
  fclose(fout);
}

// write to Gmsh format
void writeGmsh(char *meshFileName,
	       double *Node,     int N_node,
	       int    *Element,   
	       int    *ElementMark, int N_element,  
	       Face   TriangleList, int N_triangle,
	       int    verbos)           // boundary dof
{
  char meshFile[60];
  FILE *fout;
  int i,j, idx;
  double x,y,z;
  
  sprintf(meshFile,"%s.msh",meshFileName);
  if(verbos)
    printf("Writing the mesh to FILE %s with Gmsh Format.\n", meshFile);

  fout = fopen(meshFile, "w");

  // 0. write the mesh file head
  fprintf(fout, "$MeshFormat\n"); 
  fprintf(fout, "2 0 8\n");
  fprintf(fout, "$EndMeshFormat\n" );

  if(Node == NULL)
    {
      printf("The Mesh is empty!\n");
    }
  else
    {
      // 1. write the nodes' coordinates
      fprintf(fout, "$Nodes\n");
      fprintf(fout, "%d\n", N_node);
      for(i = 0; i < N_node; i++)
	{
	  idx = 3*i;
	  fprintf(fout, "%d %lf %lf %lf \n", 
		  i+1, Node[idx], Node[idx+1], Node[idx+2]);
	}
      fprintf(fout, "$EndNodes\n");
      fprintf(fout, "\n");   // empty line

      // 2. write the all the elements, which includeing the boundary face
      fprintf(fout, "$Elements\n"); 
      fprintf(fout, "%d\n", N_triangle + N_element);

      // 2.1 first boundary surface, optional
      int count = 0;
      Face *the_face = TriangleList._next;   // the head node is empty!
      while(the_face != NULL){
	fprintf(fout, "%d 2 1 %d", the_face->_id+1, the_face->_boundary_mark);
	fprintf(fout," %d %d %d\n",    the_face->_vtx[0]+1,
		the_face->_vtx[1]+1, the_face->_vtx[2]+1);

	  the_face = the_face->_next;  // go to the next triangle
	count ++;
      }
      fprintf(fout, "\n");   // empty line

      if(count != N_triangle) 
	printf("Error in writing the triangle data into gmsh format!\n");

      // 2.2 then the tetrahedron
      for(i = 0, idx = 0; i < N_element; i++, idx += 4)
	{
	  fprintf(fout, "%d 4 1 %d", N_triangle + i + 1, ElementMark[i]);

	  fprintf(fout, " %d %d %d %d \n", 
		  Element[idx + 0] + 1, 
		  Element[idx + 1] + 1, 
		  Element[idx + 2] + 1, 
		  Element[idx + 3] + 1);
	}


      fprintf(fout, "$EndElements\n");

      if(verbos)
	printf("%d Nodes, %d Triangles and %d Tetrahedron have been write into mesh file %s.\n", N_node, N_triangle, N_element, meshFileName);

    }
  // 3.end of the output
  fclose(fout);
}


// load mesh from mesh data, which is stored as in AFEPack format
void readAFEPackMesh(char *meshFileName,
		     double **Node,    int *N_node,     //
		     int    **Element, int *N_element,
		     int verbose)
{
  FILE *fin;
  int i,idx;
  int j,k,l,m,n;
  int dumid,dummy;
  double x,y,z;
  
  fin = fopen(meshFileName, "r");
  if(fin == NULL)
    {
      printf("Error in openning the mesh file: %s.\n", meshFileName);
      exit(1);
    }

  // 1. read nodes' coordinates
  fscanf(fin, "%d", N_node);
  *Node =(double *)malloc((*N_node)*3*sizeof(double));
  if (*Node == NULL) // fail to malloc memory spaces
    {
      printf("Fail to reserve memory for Node.\n");
      exit(1);
    }
  for(i = 0, idx = 0; i < *N_node; i++, idx += 3)
    {
      fscanf(fin, "%lf %lf %lf", &x, &y, &z);
      (*Node)[idx+0] = x;
      (*Node)[idx+1] = y;
      (*Node)[idx+2] = z;
    }

  // skip the 0-dim geometry object in AFEPack mesh format
  fscanf(fin, "%d", &dummy);
  for(i = 0; i < *N_node; i++)
    {
      fscanf(fin, "%d", &dummy);
      fscanf(fin, "%d %d", &dummy, &dummy);
      fscanf(fin, "%d %d", &dummy, &dummy);
      fscanf(fin, "%d", &dummy);
    }

  // 2. read the edges
  fscanf(fin, "%d", &l);
  //  *Edge            =(int *)malloc((*N_edge)*2*sizeof(int));  
  for(i = 0, idx = 0; i < l; i++, idx += 2)
    {
      fscanf(fin, "%d", &dummy);
      fscanf(fin, "%d %d %d", &dummy, &j, &k);
      fscanf(fin, "%d %d %d", &dummy, &dummy, &dummy);
      fscanf(fin, "%d", &dummy);
      //      (*Edge)[idx+0] = j;      // AFEPack numbered geometry from 0
      //      (*Edge)[idx+1] = k;      // the same with memory mesh model
    }

  // skip the triangles
  fscanf(fin, "%d", &dumid);
  for(i = 0; i < dumid; i++)
    {
      fscanf(fin, "%d", &dummy);
      fscanf(fin, "%d %d %d %d", &dummy, &j, &k, &l);
      fscanf(fin, "%d %d %d %d", &dummy, &j, &k, &l);
      fscanf(fin, "%d", &dummy);
    }

  // 3. read the elements
  fscanf(fin, "%d", N_element);
  *Element         =(int *)malloc((*N_element)*4*sizeof(int));  
  if (*Element == NULL) // fail to malloc memory spaces
    {
      printf("Fail to reserve memory for Element.\n");
      exit(1);
    }
  for(i = 0, idx = 0; i < *N_element; i++, idx += 4)
    {
      fscanf(fin, "%d", &dummy);
      fscanf(fin, "%d %d %d %d %d", &dummy, &j, &k, &m, &n);
      fscanf(fin, "%d %d %d %d %d", &dummy, &l, &l, &l, &l);
      fscanf(fin, "%d", &dummy);
      (*Element)[idx+0] = j;
      (*Element)[idx+1] = k;
      (*Element)[idx+2] = m;
      (*Element)[idx+3] = n;
    }

  if(verbose)
    printf("Read %d Nodes and %d Elements from %s in AFEPack format!\n",
	   *N_node, *N_element, meshFileName);

}


/**
 * Calculate the boundary Triangles
 * 
 */
void build_mesh_boundary(int **Triangle, int **TriangleEdge, 
 			 int **TriangleElement, int *N_triangle,
			 int *Element, int N_element,
			 int *Edge,    int N_edge,
			 int verbose)
{
  int elem, tri, eg, idx, iFace,j, pos;
  int Face_end[4][3];
  Face_end[0][0] = 1;  Face_end[0][1] = 2;  Face_end[0][2] = 3;
  Face_end[1][0] = 0;  Face_end[1][1] = 2;  Face_end[1][2] = 3;
  Face_end[2][0] = 0;  Face_end[2][1] = 1;  Face_end[2][2] = 3;
  Face_end[3][0] = 0;  Face_end[3][1] = 1;  Face_end[3][2] = 2;
  
  if(verbose)
    printf("Building the Face informations ... \n");

  int N_face = 0;
  int *Face = (int *) malloc(3*2*N_edge*sizeof(int));
  int *ElementFace = (int *) malloc(4*N_element*sizeof(int));
  for(j = 0; j < 4*N_element; j++) ElementFace[j] = 0;

  // step 1: book keep the element face
  int v0, v1, v2;
  int u0, u1, u2;
  for(elem = 0; elem < N_element; elem ++){
    for(iFace = 0; iFace < 4; iFace++){

      // vertrice of face iFace
      v0 = Element[4*elem + Face_end[iFace][0]];
      v1 = Element[4*elem + Face_end[iFace][1]];
      v2 = Element[4*elem + Face_end[iFace][2]];

      for(idx = 0; idx < N_face; idx++){
	u0 = Face[3*idx+0];
	u1 = Face[3*idx+1];
	u2 = Face[3*idx+2];

	if(((u0 == v0) && (u1 == v1) && (u2 == v2)) ||
	   ((u0 == v0) && (u1 == v2) && (u2 == v1)) ||
	   ((u0 == v1) && (u1 == v0) && (u2 == v2)) ||
	   ((u0 == v1) && (u1 == v2) && (u2 == v0)) ||
	   ((u0 == v2) && (u1 == v0) && (u2 == v1)) ||
	   ((u0 == v2) && (u1 == v1) && (u2 == v0))){ 
	  break;   // find it, do nothing but exit with idx
	} 

      }

      ElementFace[4*elem + iFace] = idx;  // add face information

      if(idx == N_face){    // when the face(v0,v1,v2) do not exist
	Face[3*idx+0] = v0; // a new face is needed 
	Face[3*idx+1] = v1; 
	Face[3*idx+2] = v2;
	N_face ++;   
      }


    }  // finish all 4 faces in element
  }    // finish all element

  // Now, there are N_face in the list
  // step 2: we can then fill the FaceElement informations

  if(verbose)
    printf("Building the Face-Element connections ... \n");

  int *FaceElement = (int *) malloc(2*N_face*sizeof(int));
  int *FaceElemIdx = (int *) malloc(N_face*sizeof(int));
  for(j = 0; j < N_face; j++){
    FaceElemIdx[j] = 0;
    FaceElement[2*j+0] = 0;
    FaceElement[2*j+1] = 0;
  }

  int face_gdx;
  for(elem = 0; elem < N_element; elem++){
    for(iFace = 0; iFace < 4; iFace++){

      face_gdx = ElementFace[4*elem + iFace];
      
      pos = 2*face_gdx + FaceElemIdx[face_gdx];    
      FaceElement[pos] = elem;   // elem is now a neighbour of face face_gdx
      
      FaceElemIdx[face_gdx]++;  // point to the next position
      if(FaceElemIdx[face_gdx] > 2){
	printf("Wrong information with Index of %d is %d.\n", 
	       face_gdx, FaceElemIdx[face_gdx]);
      }
      
    }
  }

  // step 3: count the boundary triangle with FaceElement[*][1] == 0
  *N_triangle = 0;
  for(j = 0; j < N_face; j++){
    if(FaceElement[2*j+1] == 0) {
      (*N_triangle)++;
    }
  }

  if(verbose)
    printf("Retriving the %d boundary triangles ... \n", *N_triangle);

  // step 4: calculate the boundary triangles
  (*Triangle)     = malloc(3*(*N_triangle)*sizeof(int));
  (*TriangleEdge) = malloc(3*(*N_triangle)*sizeof(int));
  (*TriangleElement) = malloc((*N_triangle)*sizeof(int));
  for(j = 0, idx = 0; j < N_face; j++){
    if(FaceElement[2*j+1] != 0) {
      continue; // skip non-boundary faces
    }

    // add three nodes to triangle idx
    (*Triangle)[3*idx+0] = Face[3*j+0];
    (*Triangle)[3*idx+1] = Face[3*j+1];
    (*Triangle)[3*idx+2] = Face[3*j+2];

    (*TriangleElement)[idx] = FaceElement[2*j+0]; //pick up the left neighbour
    
    // find the global index of the edge is tedious
    for(eg = 0; eg < 3; eg++){
      v0 = Face[3*j + (eg+1)%3];
      v1 = Face[3*j + (eg+2)%3];
      // need to find the global position of edge[v0 v1]
      for(pos = 0; pos < N_edge; pos++){
	u0 = Edge[2*pos + 0];
	u1 = Edge[2*pos + 1];
	if(((u0 == v0) && (u1 == v1)) ||
	   ((u0 == v1) && (u1 == v0))){
	  break;    // must be found
	}
      }
      if(pos == N_edge){
	printf("Something wrong!\n");
	exit(1);
      }

      (*TriangleEdge)[3*idx + eg] = pos;
    }
	
    idx ++;
    
  }

  if(Face != NULL) free(Face);
  if(FaceElement != NULL) free(FaceElement);
  if(FaceElemIdx != NULL) free(FaceElemIdx);

}



/**
 * calculate the Edges' index for each element
 * based on the arrays Element and Edge
 */

void build_element_edge(int *ElementEdge,  // output
			int *Element, int N_element,
			int *Edge, int N_edge,
			int verbose)
{
  int elem, eg, j;
  int v_start, v_end, v1, v2;
  int idx;

  int Edge_end[6][2];                   // in local edge index
  Edge_end[0][0]=0; Edge_end[0][1]=1;
  Edge_end[1][0]=0; Edge_end[1][1]=2;
  Edge_end[2][0]=0; Edge_end[2][1]=3;
  Edge_end[3][0]=1; Edge_end[3][1]=2;
  Edge_end[4][0]=1; Edge_end[4][1]=3;   // maybe should be reversed ||
  Edge_end[5][0]=2; Edge_end[5][1]=3;
  
  if(verbose)
    printf("Calculating the ElementEdge ... \n");
  for( elem = 0; elem < N_element; elem ++){
    for( eg = 0; eg < 6; eg ++){

      idx = 4*elem + Edge_end[eg][0];
      v_start = Element[idx];
      idx = 4*elem + Edge_end[eg][1];
      v_end   = Element[idx];

      // go to find the edge[v_start v_end]'s global position
      for(idx = 0; idx < N_edge; idx++){
	v1 = Edge[2*idx+0];
	v2 = Edge[2*idx+1];
	if(((v_start == v1) && (v_end == v2)) ||
	   ((v_start == v2) && (v_end == v1))){ // found
	  ElementEdge[6*elem+eg] = idx;
	  break;   // if found, go to the next edge of the element
	}
      }  // finish search Edge list

    }  // finish 6 edges of single element
  }   // finish of all element

}


/**
 * 读取Gmsh 网格文件信息
 * 
 */
int readGmshHead(char *meshFileName,
		 int  *N_node,
		 int  *N_triangle, int *n_flag_per_tri,
		 int  *N_element, int *n_flag_per_elem,
		 int   verbose)
{
  char meshFile[60], dummy[60];
  double version;
  int data_type, sizeof_double;

  int inode=0, iedge=0, itri=0, ielem=0;
  int nnode=0, nedge=0, ntri=0, nelem=0;

  int igeom = 0, ngeom = 0, nvtx = 0;
  int id, type, nflag;  // fist infomation for any element
  int i, j, k, flag[4], vtx[4]; // reserve 4 space for flag and nodes

  double x,y,z;
  FILE *f_in;
  sprintf(meshFile,"%s.msh",meshFileName);
  f_in = fopen(meshFile, "r");
  if(f_in == NULL){
    printf("The mesh file %s do not exist!\n", meshFile);
    return 1;
  }
    

  //skip the first line
  fscanf(f_in, "%s", dummy);
  fscanf(f_in, "%lf %d %d", &version, &data_type, &sizeof_double); 
  fscanf(f_in, "%s", dummy);

  if(verbose)
    printf("Gmsh: format version %f with flag %d and size of double being %d.\n", 
	   version, data_type, sizeof_double);

  // read Node information
  fscanf(f_in, "%s", dummy);
  fscanf(f_in, "%d", &nnode);
  for (inode = 0; inode < nnode; inode++)
    fscanf(f_in, "%d %lf %lf %lf", &id, &x, &y, &z);  // ingor the data
  fscanf(f_in, "%s", dummy);  // skip $EndNodes

  // read all kinds of elements:
  fscanf(f_in, "%s", dummy);   // skip $Elements
  fscanf(f_in, "%d", &ngeom);  
  for (igeom = 0; igeom < ngeom; igeom++){
    fscanf(f_in, "%d %d %d", &id, &type, &nflag);  // ingor the data
    switch(type){
    case 15:{  // node
      nvtx = 1;  // skip it
      break;
    }
    case 1:{   // edge
      nedge ++;
      nvtx = 2;
      break;
    }
    case 2:{   // triangle
      ntri ++;
      nvtx = 3;
      *n_flag_per_tri = nflag;
      break;
    }
    case 4:{    // tethedron
      nelem ++;
      nvtx = 4;
      *n_flag_per_elem = nflag;
      break;
    }

    default:
      printf("Sorry, we can only treat with tetrahedron mesh currently.\n");
      
    }
    for(i = 0; i < nflag; i++) fscanf(f_in, "%d", flag + 0);
    for(i = 0; i < nvtx;  i++) fscanf(f_in, "%d", vtx  + i);
  }
  fscanf(f_in, "%s", dummy);    // skip $EndElements

  fclose(f_in);

  // save the information for the number of nodes
  *N_node = nnode;             
  //  *N_edge = nedge;
  *N_triangle = ntri;
  *N_element = nelem;

  if(verbose)
    printf("There are %d nodes, %d edges, %d boundary tiangles and %d elements.\n", 
	   *N_node, nedge, *N_triangle, *N_element);
  
  return 0;

}

/**
 * read the concrete data
 * 
 */

int readGmshData(char *meshFileName,
		 double  *Node, int N_node,
		 int *Triangle, int N_triangle,
		 int *TriangleFlag, int n_flag_per_tri,
		 int  *Element, int N_element,
		 int *ElementFlag,  int n_flag_per_elem,
		 int   verbose)
{
  char meshFile[60], dummy[60];
  double version;
  int data_type, sizeof_double;

  int inode=0, iedge=0, itri=0, ielem=0;
  int nnode=0, nedge=0, ntri=0, nelem=0;

  int igeom = 0, ngeom = 0, nvtx = 0;
  int id, type, nflag;  // fist infomation for any element
  int i, j, k, flag[4], vtx[4]; // reserve 4 space for flag and nodes

  double x,y,z;
  FILE *f_in;
  sprintf(meshFile,"%s.msh",meshFileName);
  f_in = fopen(meshFile, "r");
  if(f_in == NULL){
    printf("The mesh file %s do not exist!\n", meshFile);
    return 1;
  }
    

  //skip the first line
  fscanf(f_in, "%s", dummy);
  fscanf(f_in, "%lf %d %d", &version, &data_type, &sizeof_double); 
  fscanf(f_in, "%s", dummy);

  // read Node information
  fscanf(f_in, "%s", dummy);
  fscanf(f_in, "%d", &nnode);
  if(nnode != N_node){
    printf("The wrong mesh informations.\n");
    return 1;  //fail
  }
  // otherwise, keep reading
  for (inode = 0; inode < nnode; inode++){
    fscanf(f_in, "%d %lf %lf %lf", &id, &x, &y, &z);  // ingor the data
    Node[3*inode + 0] = x;
    Node[3*inode + 1] = y;
    Node[3*inode + 2] = z;
  }
  fscanf(f_in, "%s", dummy);  // skip $EndNodes

  // read all kinds of elements:
  fscanf(f_in, "%s", dummy);   // skip $Elements
  fscanf(f_in, "%d", &ngeom);
  nedge = 0; ntri = 0; nelem = 0;
  for (igeom = 0; igeom < ngeom; igeom++){
    fscanf(f_in, "%d %d %d", &id, &type, &nflag);  // element information
    // ignore the flags
    for(i = 0; i < nflag; i++) fscanf(f_in, "%d", &flag[i]);
    // get the vertices 
    switch(type){    
    case 15:{  // node
      nvtx = 1;
      fscanf(f_in, "%d", &flag[0]); // skip the vertex;
      break;
    }
    case 1:{   // edge
      fscanf(f_in, "%d %d", &vtx[0], &vtx[1]);
      //     Edge[2*nedge + 0] = vtx[0] - 1;
      //      Edge[2*nedge + 1] = vtx[1] - 1;
      nedge ++;
      break;
    }
    case 2:{   // triangle
      for(i = 0; i < nflag; i++){
	TriangleFlag[ntri*n_flag_per_tri + i] = flag[i];
      }
      fscanf(f_in, "%d %d %d", &vtx[0], &vtx[1], &vtx[2]);
      Triangle[3*ntri + 0] = vtx[0] - 1; 
      Triangle[3*ntri + 1] = vtx[1] - 1;
      Triangle[3*ntri + 2] = vtx[2] - 1;
      ntri ++;
      break;
    }
    case 4:{    // tethedron
      for(i = 0; i < nflag; i++){
	ElementFlag[nelem*n_flag_per_elem + i] = flag[i];
      }
      fscanf(f_in, "%d %d %d %d", vtx + 0, vtx + 1, vtx + 2, vtx + 3);
      Element[4*nelem + 0] = vtx[0] - 1;
      Element[4*nelem + 1] = vtx[1] - 1;
      Element[4*nelem + 2] = vtx[2] - 1;
      Element[4*nelem + 3] = vtx[3] - 1;
      nelem ++;
      break;
    }

    default:
      printf("Sorry, type=%d on geometry %d,we can only treat with tetrahedron data currently.\n", type, igeom);
      
    }
  }
  fscanf(f_in, "%s", dummy);    // skip $EndElements

  fclose(f_in);

  if((ntri != N_triangle) || (nelem != N_element)){
    printf("!!! Warning: mesh data may be incorrect. !!!\n");
  }

  if(verbose)
    printf("Mesh loaded from %s.\n", meshFile);

  return 0;  // not fail

}		  


/**
 * 根据指定的边界标记保存网格到gmsh格式
 * 
 */
void write_gmsh_marked(char *mesh_prefix, 
		      double *Node, int N_node,
		      int *Element, int N_element,
		      int *ElementFlag, int n_flag_per_elem, 
		      int flag_idx)  // according to flag[idx]
{
  /**
   * 先统计单元边界标记的种类和数量
   * 
   */
  int n_bnd_mark = 0;
  int bnd_mark[20];  // reserve more spaces
  int cnt_mark[20];  // howmany element belong to specific mark
  int i,j,k,pos,flag;

  //  int n_flag_per_elem = N_element_flag/N_element;

  // initialize
  for(i = 0; i < 20; i++){
    bnd_mark[i] = 0;
    cnt_mark[i] = 0;
  }

  for(i = 0; i < N_element; i++){
    flag = ElementFlag[n_flag_per_elem*i + flag_idx];
    for(pos = 0; pos < n_bnd_mark; pos++){
      if(flag == bnd_mark[pos]){ 
	cnt_mark[pos] ++; // add one more element to this mark
	break;
      }
    }
    if(pos == n_bnd_mark){   // new boundary mark
      bnd_mark[n_bnd_mark++] = flag;
      cnt_mark[pos] = 1;  // It must be the first element in this mesh
    }
  }
  printf("There are %d different mark for the element.\n", n_bnd_mark);

  //
  char fileName[60];
  int *submesh, *subflag;
  int elem_pos, nn=n_flag_per_elem;
  for(i = 0; i < n_bnd_mark; i++){

    // 1. allocate space first
    submesh = (int *) malloc(4 *cnt_mark[i]*sizeof(int));
    subflag = (int *) malloc(nn*sizeof(int));

    // 2. fill in element data
    for(j = 0, elem_pos = 0; j < N_element; j++){
      flag = ElementFlag[n_flag_per_elem*j + flag_idx];
      if(flag != bnd_mark[i]) continue; //skip other element
      // copy the data
      submesh[4*elem_pos + 0] = Element[4*j + 0];
      submesh[4*elem_pos + 1] = Element[4*j + 1];
      submesh[4*elem_pos + 2] = Element[4*j + 2];
      submesh[4*elem_pos + 3] = Element[4*j + 3];
	subflag[elem_pos] = ElementFlag[nn*j + 0];      
      elem_pos ++;
    } 

    // 3. write it out with all the nodes 
    sprintf(fileName, "%s_%d", mesh_prefix, bnd_mark[i]);
    Face TriangleList = {-1, -1,-1,-1, -1,-1,-1, -1,-1, -1, NULL};
    writeGmsh(fileName, 
	      Node, N_node, submesh, subflag, elem_pos,
	      TriangleList, 0, 1);

    // 4. free the space for the current submesh
    if(submesh != NULL) free(submesh);
    if(subflag != NULL) free(subflag);

  }   // finish treating mark[i]

}

/**
 * 计算单纯形的体积
 * Assume in 3D space
 */

double volume_tet(double *v0,
		  double *v1,
		  double *v2,
		  double *v3){
  return ((v1[0] - v0[0])*(v2[1] - v0[1])*(v3[2] - v0[2])
	  + (v1[1] - v0[1])*(v2[2] - v0[2])*(v3[0] - v0[0])
	  + (v1[2] - v0[2])*(v2[0] - v0[0])*(v3[1] - v0[1])
	  - (v1[0] - v0[0])*(v2[2] - v0[2])*(v3[1] - v0[1])
	  - (v1[1] - v0[1])*(v2[0] - v0[0])*(v3[2] - v0[2])
	  - (v1[2] - v0[2])*(v2[1] - v0[1])*(v3[0] - v0[0]))/6.;  
}


/**
 * 修改节点的排列顺序使得所有四面体得到的Jacobian值为正
 * 
 */

void reorder_node(double *Node, int *Element, int N_element){
  int i, l, idx, count = 0;
  for(i = 0; i < N_element; i++){
    double area = volume_tet(&Node[3*Element[4*i + 0]],
			     &Node[3*Element[4*i + 1]],
			     &Node[3*Element[4*i + 2]],
			     &Node[3*Element[4*i + 3]]);
    if(area < 0){
      l = Element[4*i + 0];
      Element[4*i + 0] = Element[4*i + 1];
      Element[4*i + 1] = l;
      count ++;
    }   
  }
  printf("Converted %d elements.\n", count);
}


/**
 * build edge info from simplest mesh,
 * which only the nodes' position and the element vertices are known.
 * Assuming that the memory for output value have been allocated
 */
void build_edge(int *Edge, int *N_edge, int *ElementEdge,   // output
		int *Element, int N_element, int N_node)     // input mesh
{
  int i_elem, i_eg, n_eg, i,j, k,l;
  int v_start, v_end, v1, v2;
  int idx, pos, eg_pos;

  //  if(verbose) printf("Calculating the ElementEdge ... \n");

  /**
   * This is a optimized methodology for tetrahedron mesh:
   * if mat(i,j) == nnz, then node i and node j are connected
   *
   */
  int nnz;
  int *row_ptr, *col_idx;
  sparse_pattern_nodal(&row_ptr, &col_idx, &nnz,
		       Element, N_node, N_element);

  // 1. generate edge indices
  int *EdgeIdx = (int *) malloc(nnz*sizeof(int));
  for(i = 0; i < nnz; i++) EdgeIdx[i] = 0;  // initialize

  for(i = 0, eg_pos = 0; i < N_node; i++){
    for(pos = row_ptr[i]; pos < row_ptr[i+1]; pos++){
      j = col_idx[pos];
      if(i == j){  // diagnal element, only add once
	add_sparse_int(row_ptr, col_idx, EdgeIdx,
		       i, i, -1); // the diagnal is invalid
      }
      else if(i < j){ // it must be symmetric
	add_sparse_int(row_ptr, col_idx, EdgeIdx,
		       i, j, eg_pos);
	add_sparse_int(row_ptr, col_idx, EdgeIdx,
		       j, i, eg_pos); // add this for search table
	eg_pos ++;
      }
    }
  }

  // make sure the connectivity matrix is symmetric
  test_symm(row_ptr,col_idx,EdgeIdx,N_node);

  // 2. Fill in for Edge
  *N_edge = eg_pos;  
  for(i = 0; i < N_node; i++){
    for(pos = row_ptr[i]; pos < row_ptr[i+1]; pos++){
      j = col_idx[pos];
      if(i >= j) continue;    // only use the upper triangle
      Edge[2*EdgeIdx[pos] + 0] = i;   // the first node
      Edge[2*EdgeIdx[pos] + 1] = j;   // is smaller
    }
  }

  // 3. Fill in for EdgeElement
  for(i_elem = 0; i_elem < N_element; i_elem ++){
    for(i_eg = 0; i_eg < 6; i_eg ++){
      // vertex of eg i_eg in element i_element
      i = Element[4*i_elem + Edge_end[i_eg][0]];
      j = Element[4*i_elem + Edge_end[i_eg][1]];
      for(pos = row_ptr[i]; pos < row_ptr[i+1]; pos++){
	if(col_idx[pos] == j) break;
      }

      // check the validaty, the edge must be found in row i
      if(pos == row_ptr[i+1]){
	printf("Something wrong when building the edge!\n");
      }

      // save the edge informations
      ElementEdge[6*i_elem + i_eg] = EdgeIdx[pos];

    }  // finish 6 edges
  }   // finish all elements
  
    // free sparse matrix memory
  if(row_ptr != NULL) free(row_ptr);   
  if(col_idx != NULL) free(col_idx);
  if(EdgeIdx != NULL) free(EdgeIdx);

}


/**
 * Build face informations
 * only the interface and boundary triangle 
 * need to be returned 
 */
// use insert algorithm and linke list
// which keep increasing index
// It could reduce the search time siginificantly
// the face is sorted according to three vertices


//
// 2. search the_face with verticex(i,j,k) in FaceList, i = min(i,j,k); 
// if found, return the pointer in the list
// else, return the value of *row_ptr[i+1]
//
Face * find_face(Face *FaceList, int i, int j, int k){
  Face *the_face;
  int v0,v1,v2; // the vertices number in the list 
  
  // only search in a limited region
  the_face = FaceList[i]._next;  // skip the head nodes
  for( ; the_face != FaceList + (i+1); the_face = the_face->_next){

    v0 = the_face->_vtx[0];  
    if(v0 != i) {   // it must be i!
      printf("There must be error in sorting the Face list!\n");
    }

    v1 = the_face->_vtx[1];
    v2 = the_face->_vtx[2];
    
    if( ((v1 == j) && (v2 == k)) || ((v1 == k) && (v2 == j)) ) { 
      break;  // find it, do nothing
    }
       
  }

  // case 1: find (the_face != FaceList+i+1), 
  //   then you can get further infomation from it
  // case 2: not find (the_face == FaceList+i+1)
  // it is easy to add a new face at FaceList+i
  return the_face;    
}

//
// 3. the main function to build all triangle faces
//
void build_face(Face *FaceList,  int *N_face,   //output, insert all faces
		int *Element, int *ElementEdge, int N_element, // input
		int *Edge,    int N_edge, int N_node)  //input
{
  int i,j,k,l;
  int vtx[3],ldx[3];
  int pos,idx;
  int i_elem, i_face, i_node, n_face = 0;
  /**
   * 0. initialize the FaceList's head nodes
   * each node is associated with a list of faces
   */

  for(i_node = 0; i_node < N_node; i_node++){  // initialize
    FaceList[i_node]._id = -1;

    FaceList[i_node]._vtx[0] = -1;
    FaceList[i_node]._vtx[1] = -1;
    FaceList[i_node]._vtx[2] = -1;

    FaceList[i_node]._bnd[0] = -1;
    FaceList[i_node]._bnd[1] = -1;
    FaceList[i_node]._bnd[2] = -1;

    FaceList[i_node]._nei[0] = -1;
    FaceList[i_node]._nei[1] = -1;

    FaceList[i_node]._boundary_mark = -1;

    FaceList[i_node]._next = &FaceList[i_node + 1]; // one more space due to this!
  }

  /**
   * 1. build the face list, sorted lexigraphically 
   * 
   */

  n_face = 0;
  int ldx_min = 0, vtx_min = 0;
  int v0,v1,v2, bnd[3];
  Face *the_face;   
  for(i_elem = 0; i_elem < N_element; i_elem++){
    for(i_face = 0; i_face < 4; i_face ++){
      // A: retrive current faces' node and vertices only
      vtx[0] = Element[4*i_elem + Face_end[i_face][0]];
      vtx[1] = Element[4*i_elem + Face_end[i_face][1]];
      vtx[2] = Element[4*i_elem + Face_end[i_face][2]];

      bnd[0] = ElementEdge[6*i_elem + Face_edge[i_face][0]];
      bnd[1] = ElementEdge[6*i_elem + Face_edge[i_face][1]];
      bnd[2] = ElementEdge[6*i_elem + Face_edge[i_face][2]];

      // B: find out the local index inside the triangle face
      
      ldx[0] = 0; ldx[1] = 1; ldx[2] = 2; 
      vtx_min = vtx[0]; ldx_min = 0;   // init the minum
      if(vtx[1] < vtx_min) {
	vtx_min = vtx[1]; ldx_min = 1;
      }
      if(vtx[2] < vtx_min) {
	vtx_min = vtx[2]; ldx_min = 2;
      }

      ldx[0] = ldx_min;          v0 = vtx[ldx[0]]; 
      ldx[1] = (ldx[0] + 1) % 3; v1 = vtx[ldx[1]];
      ldx[2] = (ldx[1] + 1) % 3; v2 = vtx[ldx[2]];

      // The index turn of edge is different from nodes'
      if(ldx[0] == 1){  // node 1 is smallest, from 0-1-2 to 2-0-1
	ldx_min = bnd[0];
	bnd[0] = bnd[2];
	bnd[2] = bnd[1]; 
	bnd[1] = ldx_min;
      }
      if(ldx[0] == 2){ // node 2 is smallest, from 0-1-2 to 1-2-0
	ldx_min = bnd[0];
	bnd[0] = bnd[1];
	bnd[1] = bnd[2]; 
	bnd[2] = ldx_min;
      }
      
      the_face = find_face(FaceList, v0, v1, v2);

      // C: postprocess after search
      if(the_face == FaceList + (v0 + 1)){ // do not find
	// generate a new face at row_ptr[v0]
	the_face = FaceList + v0; // the insert position
	Face *new_face = (Face *) malloc(sizeof(Face));
	new_face->_id = n_face++;
	new_face->_vtx[0] = v0;
	new_face->_vtx[1] = v1;
	new_face->_vtx[2] = v2;

	new_face->_bnd[0] = bnd[0];
	new_face->_bnd[1] = bnd[1];
	new_face->_bnd[2] = bnd[2];

	new_face->_nei[0] = i_elem;
	new_face->_nei[1] = -1;  // NULL currently
	
	new_face->_next = FaceList[v0]._next;   // connect to the next one in node i's list

	FaceList[v0]._next = new_face;          // and do insert

      }
      else{  //find it
	the_face->_nei[1] = i_elem;            // just add the other neighbour information
      }

      
    } // finish all 4 faces in tet
  }   // run over all tets

  *N_face = n_face;


} // end of function build_face


/**
 * retrive interface and boundary triangle from all face list
 * 
 */

void build_triangle(Face *TriangleList, int *N_triangle, //output
		    Face *FaceList, int *ElementMark, int N_node)   //input
{ 
  int i_node;
  int flag_left, flag_right, nei_right, nei_left;
  Face *the_face, *new_face; // this list is only for 
  int n_triangle = 0;
  for(i_node = 0; i_node < N_node; i_node ++){
    the_face = FaceList[i_node]._next;
    for(; the_face != FaceList + (i_node + 1); the_face = the_face->_next){
      nei_left = the_face->_nei[0];
      flag_left = ElementMark[nei_left];
      nei_right = the_face->_nei[1];
      if(nei_right != -1){   // neither boundary triangle
	flag_right = ElementMark[nei_right];
	if(flag_left == flag_right) continue; // nor interface triangle
      }

      // when excute here, either boundary triangle or interfaicial triangle
      new_face = (Face *) malloc(sizeof(Face));
      new_face->_id = n_triangle ++;
      new_face->_vtx[0] = the_face->_vtx[0];
      new_face->_vtx[1] = the_face->_vtx[1];
      new_face->_vtx[2] = the_face->_vtx[2];
      
      new_face->_bnd[0] = the_face->_bnd[0];
      new_face->_bnd[1] = the_face->_bnd[1];
      new_face->_bnd[2] = the_face->_bnd[2];

      new_face->_nei[0] = the_face->_nei[0];
      new_face->_nei[1] = the_face->_nei[1];    

      //case 1: boundary triangle
      if(nei_right == -1){  //this is the boundary
	new_face->_boundary_mark = BM_OUTER; // mark outer boundary 

      }
      else{        // case 2: must be on some interface
	new_face->_boundary_mark = BM_INTERFACE; // mark outer boundary with 
      }

      // add the new_face into TraingleList from beginning
      new_face->_next = TriangleList->_next;  
      TriangleList->_next = new_face;

    } // finish the face related with node i_node
  }  // finish all faces

  *N_triangle = n_triangle;
}

//
void mirror_mesh(double **Node, int *N_node, 
		 int **Element, int *N_elem, int axis)
{
  int i,j,e,v,pos,N_bndNode;
  int N_node0, N_elem0, N_node1, N_elem1;
  double *Node0, *Node1;
  int *NodeFlag, *bndNode, *NodeMap;
  int *Element0, *Element1;
  double x[3];

  // the number of element must doubled
  N_elem0 = *N_elem;
  Element0 = *Element;
  N_elem1 = 2*N_elem0;
  Element1 = (int *) malloc(N_elem1*4*sizeof(int));
  
  // the number of new nods should be calculated
  N_node0 = *N_node;
  Node0 = *Node;
  NodeFlag = (int *) malloc(N_node0*sizeof(int));
  for(i = 0, N_bndNode = 0; i < N_node0; i++){
    if(fabs(Node0[3*i + axis] - 0.0) > EPSDIST)
      NodeFlag[i] = 0;
    else{
      NodeFlag[i] = 1;
      N_bndNode ++;
    }
  }
  N_node1 = 2*N_node0 - N_bndNode;
  Node1 = (double *)malloc(N_node1*3*sizeof(double));

  // parepare the node map, only need to check this boundary node
  bndNode = (int *) malloc(N_bndNode*sizeof(int));
  for(i = 0, e = 0; i < N_node0; i++){
    if(NodeFlag[i] == 0) continue;
    bndNode[e++] = i;
  }
  if(N_bndNode != e) printf("wrong in building the boundary Node map");
  

  /**
   * begin to treat the new mesh
   * 
   */
  // 1. copy the old mesh
  for(i = 0; i < N_node0; i++){
    Node1[3*i + 0] = Node0[3*i + 0];
    Node1[3*i + 1] = Node0[3*i + 1];
    Node1[3*i + 2] = Node0[3*i + 2];
  }
  for(e = 0; e < N_elem0; e++){
    Element1[4*e + 0] = Element0[4*e + 0];
    Element1[4*e + 1] = Element0[4*e + 1];
    Element1[4*e + 2] = Element0[4*e + 2];
    Element1[4*e + 3] = Element0[4*e + 3];
  }

  // 2. generate the new mesh coordinates and generate the node map
  NodeMap = (int *)malloc(N_node0*sizeof(int));
  for(i = 0, pos = 0; i < N_node0; i++){
    x[0] = Node0[3*i + 0];
    x[1] = Node0[3*i + 1];
    x[2] = Node0[3*i + 2];
    // do the mirror change
    x[axis] *= -1;
    
    // if the point is not at the symmetic plane
    if(fabs(x[axis] - 0.0) > EPSDIST){
      j = N_node0 + pos;  // the new index
      pos ++;
    }
    else{ // need to find the node out
      for(v = 0; v < N_bndNode; v++){
	j = bndNode[v];
	if((fabs(x[0] - Node0[3*j + 0]) < EPSDIST) &&
	   (fabs(x[1] - Node0[3*j + 1]) < EPSDIST) &&
	   (fabs(x[2] - Node0[3*j + 2]) < EPSDIST)){
	  break; // find the node j
	}
      }
      if(v == N_bndNode) 
	printf("wrong in check the duplicated nodes for Point(%lf,%lf, %lf).\n", 
	       x[0], x[1], x[2]);
    }

    // write the node into new coordinates
    NodeMap[i] = j;
    Node1[3*j + 0] = x[0];
    Node1[3*j + 1] = x[1];
    Node1[3*j + 2] = x[2];
  }

  if(N_node0 + pos != N_node1) printf("Wrong number of new mesh!\n");

  // at last, copy the new element
  for(e = 0; e < N_elem0; e++){
    j = e + N_elem0;
    Element1[4*j + 0] = NodeMap[Element0[4*e + 0]];
    Element1[4*j + 1] = NodeMap[Element0[4*e + 1]];
    Element1[4*j + 2] = NodeMap[Element0[4*e + 2]];
    Element1[4*j + 3] = NodeMap[Element0[4*e + 3]];
  }

  // set the return value
  (*Element) = Element1;
  (*N_elem) = N_elem1;
  (*Node) = Node1;
  (*N_node) = N_node1;

  // free mem
  free(NodeMap);
  free(Element0);
  free(Node0);
  free(bndNode);
  free(NodeFlag);

  printf("There are %d Nodes and %d Elements in the new mesh.\n", *N_node, *N_elem);
}



/**
 * read the mesh information 
 * #Node and #Element
 * if succeed
 *     return 1 
 * else 
 *     return 0
 *
 * Test passed on 18 June.
 *
 */

int load_mesh_info(char *fileName, int *N_node, int *N_element)
{
  char meshFile[80];
  FILE *fin;
  sprintf(meshFile, "mesh_data/%s.info", fileName);
  fin = fopen(meshFile, "r");
  if(fin == NULL){
    printf("The mesh head file \"%s\" is not exist!\n", meshFile);
    return 0;
  }
  fscanf(fin, "%d %d", N_node, N_element);
  fclose(fin);
  return 1;
}

/**
 * load the mesh data
 * 
 * if succeed
 *     return 1 
 * else 
 *     return 0
*
 * Test passed on 18 June.
 *
 */
 
int load_mesh_data(char *fileName, 
		    double *Node, int N_node, 
		    int *Element, int N_element)
{
  char meshFile[60];
  FILE *fin;
  int v,e,i,j,m,n;
  double x,y,z;

  // 1. read the nodes' position
  sprintf(meshFile,"mesh_data/%s.node",fileName);
  fin = fopen(meshFile, "r");
  if(fin == NULL){
    printf("The mesh node file \"%s\" is not exist!\n", meshFile);
    return 0;
  }
  for(v = 0; v < N_node; v++){
    fscanf(fin, "%lf %lf %lf", &x, &y, &z);
    Node[3*v + 0] = x;
    Node[3*v + 1] = y;
    Node[3*v + 2] = z;
  }
  fclose(fin);

  // 2. read element vertices
  sprintf(meshFile,"mesh_data/%s.elem",fileName);
  fin = fopen(meshFile, "r");
  if(fin == NULL){
    printf("The mesh element file \"%s\" is not exist!\n", meshFile);
    return 0;
  }
  for(e = 0; e < N_element; e++){
    fscanf(fin, "%d %d %d %d", &i, &j, &m, &n);
    Element[4*e + 0] = i - 1;
    Element[4*e + 1] = j - 1;
    Element[4*e + 2] = m - 1;
    Element[4*e + 3] = n - 1;
  }
  fclose(fin);

  return 1;
}

/**
 * only export the mesh data to DX format
 * which is composed of xxx.node      // coordinates of nodes
 *                      xxx.elem      // start from 0
 *                      xxx.info      // number of nodes and elements
 */
void writeMeshDX(char *fileName,
		 double *Node, int N_node,
		 int *Element, int N_element)
{
  int i,e,v;
  FILE *fout;
  char file[80];
  sprintf(file, "save_data/%s.info", fileName);
  fout = fopen(file, "w");
  fprintf(fout, "%d %d ", N_node, N_element);
  fclose(fout);

  sprintf(file, "save_data/%s.node", fileName);
  fout = fopen(file, "w");
  for (v = 0; v < N_node; v ++) {
    fprintf(fout, "%lf %lf %lf \n", Node[3*v + 0], Node[3*v + 1], Node[3*v + 2]);
  }
  fclose(fout);

  sprintf(file, "save_data/%s.elem", fileName);
  fout = fopen(file, "w");
  for (e = 0; e < N_element; e++) {
    fprintf(fout, "%d %d %d %d \n", 
	    Element[4*e + 0],
	    Element[4*e + 1],
	    Element[4*e + 2],
	    Element[4*e + 3]); 
  }
  fclose(fout);
  printf("Write the mesh into openDX format!\n");
}

//

double element_volume(double *v0, double *v1, double *v2, double *v3){
  double vol = 0;
  double xe[4], ye[4], ze[4], ae[4];

  xe[0] = v0[0]; xe[1] = v1[0]; xe[2] = v2[0]; xe[3] = v3[0];
  ye[0] = v0[1]; ye[1] = v1[1]; ye[2] = v2[1]; ye[3] = v3[1];
  ze[0] = v0[2]; ze[1] = v1[2]; ze[2] = v2[2]; ze[3] = v3[2];

  // calculate all coefficients for all basis
  ae[0] =  xe[1]*(ye[2]*ze[3]-ze[2]*ye[3])
    - xe[2]*(ye[1]*ze[3]-ze[1]*ye[3]) + xe[3]*(ye[1]*ze[2]-ze[1]*ye[2]);
  ae[1] = -xe[0]*(ye[2]*ze[3]-ze[2]*ye[3])
    + xe[2]*(ye[0]*ze[3]-ze[0]*ye[3]) - xe[3]*(ye[0]*ze[2]-ze[0]*ye[2]);
  ae[2] =  xe[0]*(ye[1]*ze[3]-ze[1]*ye[3])
    - xe[1]*(ye[0]*ze[3]-ze[0]*ye[3]) + xe[3]*(ye[0]*ze[1]-ze[0]*ye[1]);
  ae[3] = -xe[0]*(ye[1]*ze[2]-ze[1]*ye[2])
    + xe[1]*(ye[0]*ze[2]-ze[0]*ye[2]) - xe[2]*(ye[0]*ze[1]-ze[0]*ye[1]);

  vol = (ae[0] + ae[1] + ae[2] + ae[3])/6.0;
  return vol;
}


double volume(double *Node, int *Element, int N_element){
  double vol = 0, vole;
  int l, i, count = 0;
  double *pt[4];
  for(i = 0; i < N_element; i++){
    pt[0] = Node + 3*Element[4*i + 0];
    pt[1] = Node + 3*Element[4*i + 1];
    pt[2] = Node + 3*Element[4*i + 2];
    pt[3] = Node + 3*Element[4*i + 3];
    vole = element_volume(pt[0],pt[1],pt[2],pt[3]);
    if(vole < 0){
      l = Element[4*i + 0];
      Element[4*i + 0] = Element[4*i + 1];
      Element[4*i + 1] = l;
      count ++;
    }
    vol += fabs(vole);
  }
  return vol;
}
